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ABSTRACT 

The non-local hydrodynamic moment equations for compressible convection are compared to 
numerical simulations. Convective and radiative flux typically deviate less than 20% from the 
3D simulations, while mean thermodynamic quantities are accurate to at least 2% for the cases 
we have investigated. The moment equations are solved in minutes rather than days on standard 
workstations. We conclude that this convection model has the potential to considerably improve 
the modelling of convection zones in stellar envelopes and cores, in particular of A and F stars. 

Subject headings: hydrodynamics — turbulence — convection — stars: interiors — stars: atmospheres 



1. Introduction 

In a wide range of astrophysical problems heat 
is transported by both radiation and convection. 
Examples include stellar envelopes and stellar 
cores, where convection may be coupled with ro- 
tation, pulsation, and diffusion. To be useful, a 
convection model must be both manageable and 
reliable. Local convection models satisfy both cri- 
teria but they are restricted to regions of strong to 
moderately strong convection, while in most cases 
one needs to model moderately strong to weak 
convection. 

Numerical simulations (Sofia and Chan (1984); 
Freytag et al. (1996), and Atroshchenko and 
Gadun (1994); Kim and Chan (1998); Stein and 
Nordlund (1998)) have been applied to 2D and 
3D stellar atmospheres. Their successful repro- 
duction of spectral line profiles and solar granu- 
lation statistics has proved their reliability. How- 
ever, they are restricted to layers near the stellar 
surface (Kupka 1999a) because of the huge ther- 
mal time scales characterizing stellar interiors. In 
addition, they currently are too expensive for ap- 
plications such as non-linear stellar pulsation cal- 
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culations for RR Lyrae stars (Feuchtinger 1998) 
or for use in everyday spectrum synthesis of large 
wavelength ranges or large parameter sets for A- 
M stars. Hence, they do not satisfy the criterion 
of manageability for the whole range of problems 
where convection occurs. 

Convection models based on the non-local, hy- 
drodynamic moment equations provide a possible 
alternative. They describe the moments of an en- 
semble average of the basic fields: velocity U, tem- 
perature T, and density p (or pressure P). Dy- 
namic equations for the moments are derived di- 
rectly from the fully compressible Navier-Stokes 
equations (NSE, Canuto (1997)). The "ensemble" 
used in this averaging process consists of realiza- 
tions of solutions of the NSE (Xiong 1986, 1997; 
Canuto 1992, 1993, 1997; Canuto and Dubovikov 
1998), or of "convective elements" of different ve- 
locity and temperature (Grossman et al. 1993; 
Grossman 1996). However, the moment equations 
entail higher order moments and thus require clo- 
sure assumptions. To obtain a closed set of equa- 
tions most convection models have invoked a mix- 
ing length. The models derived in Canuto (1992, 
1993, 1997) and in Canuto and Dubovikov (1998) 
avoid a mixing length by providing a fully non- 
local set of dynamic equations for the second order 
moments. 
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Recently, Kupka (1999a) has presented the first 
comparison of a variant of these convection models 
with numerical simulations of compressible con- 
vection for a stellar-like scenario. He used the 
downgradient approximation (DGA) for the third 
order moments. To avoid its shortcomings, here 
we use a more complete model. It is based on a 
model introduced in Canuto (1992). The latter 
was successfully applied to the convective bound- 
ary layer of the terrestrial atmosphere (Canuto et 
al. 1994). We describe the physical scenario for 
our comparison and present results for two sample 
problems. We also consider the potential of such 
models for application to envelope convection in 
A and F stars. 

2. The physical scenario 

In Muthsam ct al. (1995, 1999) the fully com- 
pressible NSE were solved for a 3D plane parallel 
geometry with a constant gravity g pointing to 
the bottom of a simulation box. Periodic bound- 
ary conditions were assumed horizontally. A con- 
stant temperature T was prescribed at the top 
and a constant input flux at the bottom of the 
box (dT/dz — const., z denotes the vertical di- 
rection, i.e. top to bottom). Top and bottom were 
taken to be impenetrable and stress free. A perfect 
gas law was assumed and radiation was treated 
in the diffusion approximation. Stable and un- 
stable layers were defined through dT/dz, which 
initially was a piecewise linear function in units of 
the adiabatic gradient: dT/dz = b(z)(dT/dz) a , < i, 
where b(z) = 1 + (V - V ad )/V ad . To define the 
stability properties of the layers a Rayleigh num- 
ber Ra for one zone where b{z) — const, and > 1 
(Muthsam ct al. 1995, 1999) is specified together 
with a Prandtl number Pr and an initial temper- 
ature contrast. This yields the radiative conduc- 
tivity K ra d(z). Both Pr and if ra d are kept fixed 
and place regions of stable and unstable stratifica- 
tion at different vertical locations in the simulation 
box. A similar approach was used by Cattaneo et 
al. (1991); Chan and Sofia (1996); Hurlburt et al. 

(1994) ; Porter and Woodward (1994); Singh et al. 

(1995) and others. The comparison to simulations 
with a prescribed viscosity avoids the need to use 
a subgrid scale model. According to our numerical 
experiments with the moment equations, molecu- 
lar viscosity can decrease the efficiency of convec- 
tion as measured by F conv /F to tai by up to 15% 



and smooths out the numerical solution in stably 
stratified regions (Kupka 1999a, b). 

3. A convection model based on the hy- 
drodynamic moment equations 

Using the Reynolds stress approach, Canuto 
(1992, 1993) has derived a convection model which 
consists of four differential equations for the ba- 
sic second order moments K (turbulent kinetic 
energy), 9 2 (mean square of temperature fluctu- 
ations, i.e. thermal potential energy), w8 (F c = 
CppJ = CppwO is the convective flux), and the ver- 
tical turbulent kinetic energy w 2 . The model was 
rederived by Canuto and Dubovikov (1998, CD98) 
using a new turbulence model based on renormal- 
ization group techniques. In CD98 notation, the 
convection model reads (d t = d/dt, d z = d/dz): 

d t K + D { (K)=gaJ-e+^C u + d z (ud z K), (1) 

a t (i^) + A(^)-/?J-v 1 ^ 
+ \dz{xd*e s ) + \co, (2) 

d t J + D f (J) =(3w^+ ^ga¥ - T~gJ 

+ \d z {xd z J) + C z + ^d z {vd z J), (3) 

d t ( l -w^) + Dti^P) = -r- 1 ^ \K) 

+ 2 -gaJ - l -e+ ^C 33 + \d z (yd~^) (4) 

d t e + D f (e) = aeK^gaJ - c^K' 1 + c 3 eN 
+ d z (vd z e), N = y/j$0\, (5) 

D f (e) = d z (ew) = - l -d z \{v x + X t)d z e] (6) 

Here, a is the volume compressibility (= 1/T for 
a perfect gas), /3 is the superadiabatic gradient, 
X = Kr&d/CpP, an d v = xP r - The r's are time 
scales. We use (25a), (27b), and (28b) of CD98 
to relate the latter to the dissipation time scale 
t = 2K/e. Compressibility effects are represented 
by Cii,C0,C 3 ,C 33 given by equations (42)-(48) 



2 



of Canuto (1993). We neglected a few terms of 
Cu , Cg , C3 , C33 that were too small to contribute 
to the solution of the moment equations (Kupka 
1999b). In equation (6), v t = C^K 2 /e, and is 
a constant given by (24d) of CD98, for which we 
take the Kolmogorov constant Ko = 1.70, while 
Xt is given by the low viscosity limit of (llf ) of 
CD98. We optionally included molecular dissi- 
pation by restoring the largest (i.e. second order 
moment) terms containing the kinematic viscos- 
ity v (i.e. d z {vd z K), etc.). They are important 
when Pr is of order unity rather than zero (Kupka 
1999a,b). Hence, we included them in all the ex- 
amples shown below. For the same reason, we op- 
tionally included a term in (5) that accounts for 
molecular dissipation effects. We use c\ = 1.44, 
c 2 = 1.92, and c 3 = 0.3 where (3 < while C3 = 
elsewhere (as in CD98). The local limit of (5), 
e = if 3 / 2 /A with A = aH p and H p = P/(pg), fails 
to yield reasonable filling factors (Kupka 1999a) 
and was thus avoided. 

To calculate the mean stratification we solve 

d z (P + Pt) = -gp, (7) 
c vP d t T + pd t K = -d z (F T + F c + F k ) (8) 

where p t — pw 2 . Equations (7)-(8) are taken 
from equation (103) of Canuto (1993) (excluding 
the higher order term in his equation (104)), and 
from equation (18c) of CD98 (for the latter, c p was 
substituted to c v to account for non-Boussinesq 
effects, Canuto, priv. communication). In the sta- 
tionary limit, (8) yields F tota i = F T +F c +Fk where 
the kinetic energy flux is given by F^ — pKw and 
the radiative flux F r = —K^&dT/dz. 

Non-locality is represented by the terms Df(K), 
Df(^9 2 ), Df(J), and Df (|w 2 ), which require third 
order moments (TOMs). Using the downgradi- 
ent approach (DGA) for the TOMs Kupka (1999a) 
found qualitatively acceptable results for the con- 
vective flux and filling factors. But quantitatively 
the results were not satisfactory, which corrob- 
orated the shortcomings of the DGA found in 
Canuto et al. (1994) and in Chan and Sofia (1996). 
To improve over the DGA, one must solve the dy- 
namic equations for the TOMs (see Canuto (1992, 
1993)). We investigated both the fully time de- 
pendent case which requires 6 additional differen- 
tial equations and various approximations of their 
stationary limit which do not entail further differ- 
ential equations. Here, we consider the following 



"intermediate" model for the TOMs: we take the 
stationary limit of the dynamic equations given in 
Canuto (1993), but neglect the Boussinesq terms 
which depend on (3 (Kupka 1999b). This model 
is similar in robustness to the DGA, but yields 
significantly better results. 

For the boundary conditions wc impose T = 
const., w 2 = 0, J = 0, dK/dz = 0, de/dz = 0,¥ = 
at the top, while at the bottom dT/dz = 0, w 2 = 
0, J = 0, dK/dz = 0, de/dz = 0, dff^/dz = 0. This 
choice permits stable numerical solutions which 
are consistent with the boundary conditions for 
the numerical simulations. The equation for the 
mean pressure is constrained by varying P at the 
top such that for a given T the resulting density 
stratification is consistent with mass conservation. 

4. Results and discussion 

We have used simulation data representing two 
configurations: model 3 J, a convection zone cm- 
bedded between two strongly stable layers (Muth- 
sam et al. 1995), and model 211P, a stably strat- 
ified layer embedded between two unstable layers 
(with a small stable layer at the bottom, Muth- 
sam et al. (1999)). In both cases, radiation trans- 
ports at least 80% of the input energy and Pr = 1. 
Each unstable layer has a thickness of « 0.5-2.5 
Hp. Model 3 J has an initial adiabatic tempera- 
ture contrast of 3.5 and encompasses 4.2 H p while 
model 211P has a contrast of 6.0 and encompasses 
about 4.8 Hp. The numerical simulations were 
done for 72 x 50 x 50 grid points and have suc- 
cessfully been compared to higher resolution sim- 
ulations (125 x 100 x 100 grid points, Kupka and 
Muthsam (1999)). All simulation data shown here 
are statistical averages over many dozen sound 
crossing times. 

The moment equations were solved by centered 
finite differences on a staggered mesh. Solutions 
were calculated using 72 grid points and were 
successfully verified by comparing them with re- 
sults obtained from higher resolutions of 128 to 
512 grid points. Time integration was done by 
the Euler forward method until a stationary state 
was reached after 1.2 to 2.5 thermal time scales. 
Verification of stationarity was done by testing 
whether all errors in time derivatives were formally 
less than 10~ 8 in relative units at all grid points 
(i.e. far smaller than the truncation error) and by 
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checking the strict energy flux and mass conserva- 
tion required for stationary solutions. 

Fig. 1 compares the convective flux of model 
21 IP with two convection models: the non-local 
model with intermediate TOM and its local coun- 
terpart (the stationary, local limit of the non-local 
model, see CD98). Clearly, only the non-local 
model successfully reproduces the simulations also 
in the central overshooting region which connects 
both convection zones of model 211P. 

Fig. 2 compares the filling factor a, i.e. the rela- 
tive area in each layer covered by upwards flowing 
material, of the numerical simulations of model 
3J with the filling factor computed from the non- 
local moment equations. For the latter we used 
the prescription described in CD98. As a is an 
important topological quantity describing the in- 
homogenous nature of convection, the reasonable 
agreement found here is a very promising result. It 
illustrates the importance of improving the TOMs, 
because the DGA used in Kupka (1999a) only in- 
dicated a correct trend, while the more complete 
TOM used here also provides a closer quantitative 
agreement. The local convection model predicts a 
structureless a = 0.5. 

Fig. 3 compares the convective flux from nu- 
merical simulations for model 3J with solutions 
from the moment equations using the intermedi- 
ate TOM. Clearly, the latter has improved over the 
DGA (Kupka 1999a) both qualitatively and quan- 
titatively. Except for the overshooting region the 
DGA barely differed from the local model which 
is shown here for comparison. Though the present 
convection model provides a substantial improve- 
ment, the convective flux still falls short in the 
middle of the convection zone of 3J and the extent 
of the lower overshooting zone is underestimated. 
This may be due to neglecting effects of radia- 
tive losses on the time scales Tg and r p g as well 
as because several contributions to the TOMs are 
neglected in our "intermediate model" , or due to 
the Boussinesq treatment of the TOMs in (Canuto 
1992, 1993), or incompleteness of equation (5). 

Fig. 4 compares local length scales A = aH p as 
used in mixing length theory with the length scale 
A = K 1 - 5 /e obtained with the non-local moment 
equations for the case of model 211P. Obviously, 
there is no a that brings the length scales in agree- 
ment. Thus, there is no alternative but to solve at 
least the full equations (5)-(6) to obtain e. 



In conclusion, we have found not only quali- 
tative, but also quantitative agreement between 
numerical simulations and the non-local moment 
equations, provided one avoids the DGA for the 
TOMs and employs the stationary solution of their 
dynamic equations as suggested in appendix B of 
(Canuto 1993) and neglects Boussinesq type fac- 
tors (involving j3). Moreover, we have found the 
mean values for T, P, and p to be accurate to 
within 2% in comparison with 4% found for local 
models with optimized mixing length parameter. 
Convective and radiative flux are typically accu- 
rate to within 20%. The improvements are largest 
in regions of weak convection and in stably strat- 
ified layers. We have used the closure constants 
suggested in Canuto (1993) and CD98 except for 
en, which was increased from 0.2 to 0.5, because 
this improved the results for the planetary bound- 
ary layer (Canuto priv. communication) and en- 
hanced the numerical stability. We have not found 
the DGA to be able to yield a similar agreement 
for both 3J and 211P even when tuning the closure 
constants individually for 3 J and 2 IIP. The situ- 
ation is even worse for the local model. In Kupka 
(1999b) we will deal with this problem in detail. 
Finally, while each 3D simulation took between 
several days and several weeks on a modern work- 
station to obtain a thermally relaxed solution, the 
moment equations took a couple of minutes to half 
an hour. This holds true already for a low numer- 
ical resolution and for an explicit time integration 
method. 

The results found here are promising for the 
application of the non-local convection model in 
particular to envelope convection zones of A and 
F stars, because they feature a similar range of 
convective efficiency, thickness in terms of H p , and 
interaction among neighbouring convection zones. 
For A and F stars the thermal structure is not 
known in advance. Hence, thermal relaxation and 
thus the gain on speed in comparison with nu- 
merical simulations is essential. The computa- 
tional savings are very attractive also for prob- 
lems studied on hydrostatic time scales, when- 
ever the full information provided by a simula- 
tion is not needed. Improvements of the convec- 
tion model studied are possible, nevertheless it is a 
more promising basis for asteroseismological stud- 
ies of pulsating A and F stars (6 Set, 7 Dor, roAp, 
etc.) than local convection models and also a new 
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basis for related work such as diffusion calcula- 
tions. Detailed results for a broader range of phys- 
ical parameters, in particular for efficient convec- 
tion and deeper convection zones, and thus of high 
importance also for other types of stars, have yet 
to corroborate this study. 

I am indebted to H.J. Muthsam for permission 
to use his simulation code and data. I am grate- 
ful to V.M. Canuto and M.S. Dubovikov for dis- 
cussions on turbulent convection models. The re- 
search was performed within project P11882-PHY 
"Convection in Stars" of the Austrian Fonds zur 
Forderung der wissenschaftlichen Forschung. 
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Fig. 1. — Convective flux (model 211P) for the 
non-local convection model and for its local limit 
(mixing length A = H p ), and a numerical simula- 
tion. 

Fig. 2. - Filling factor a (model 3 J) for the non- 
local convection model and a numerical simula- 
tion. For the local convection model a = 0.5. 

Fig. 3. — Convective flux (model 3J) for the non- 
local convection model and for its local limit (mix- 
ing length A = H p ), and a numerical simulation. 

Fig. 4. Length scale A = K 1 - 5 /e (model 211P) 
from the non-local convection model compared 
with A = aH p for a = 0.7, 2.0. 
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